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The discrete Fourier transform is approximated by summing over part of the terms with corre- 
sponding weights. The approximation reduces significantly the requirement for computer memory 
storage and enhances the numerical computation efficiency with several orders without loosing ac- 
curacy. As an example, we apply the algorithm to study the three-dimensional interacting electron 
gas under the renormalized-ring-diagram approximation where the Green's function needs to be self- 
consistently solved. We present the results for the chemical potential, compressibility, free energy, 
entropy, and specific heat of the system. The ground-state energy obtained by the present calculation 
is compared with the existing results of Monte Carlo simulation and random-phase approximation. 
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I. INTRODUCTION 

For dealing with some physical problems, we need to 
take the discrete Fourier transform. Especially, many 
physical problems arc defined on lattice models. In such 
a case, we may face to the problem of Fourier transform- 
ing a function defined on the lattice to the corresponding 
reciprocal space. For most of the problems, the function 
to be transformed cannot be expressed analytically but 
given numerically. When the function needs to be given 
at a large number of discrete points within the region it 
is defined, the memory volume for storing the function 
may be too big and may even exceed the computer's stor- 
age limit. Even if the problem is within the computer's 
capability, when the transform is involved in an integral 
equation that may be solved by iterations, the function 
needs to be determined again and again in the iterations 
and the process is very time consuming. Therefore, an 
approximation scheme for the discrete Fourier transform 
that reduces the storage requirement and accelerates the 
numerical computation process without loosing the accu- 
racy is very desirable. 

The discrete Fourier transform as well as a continuous 
one is useful in solving the integral equations with con- 
volutions involved. One of the examples in the quantum 
many-body problems is to calculate the self-energy S of 
electrons p], 0], 



E(k, iu) n ) = —— ^2 v e ff(k - k', lUJ n 



iu n i)G(h! ', iu) n i) 



(1) 

where k is the momentum of electron, uj n = nT(2n + 1) 
with n an integer is the fermionic Matsubara frequency, 
T is temperature, V is the volume of the system, v e f f is 
an effective interaction between electrons, and G is the 
Green's function of electrons. In a sophisticated scheme, 
£ and G need to be determined self-consistently. After 
the Fourier transforms, in real coordinate r and imagi- 
nary time t space, Eq. ([IJ reads 



Having E(r, r) been simply calculated by Eq. ((2}, one 
then obtains iui n ) by the inverse Fourier transforms. 

Earlier works dealt with Eq. ((T|) by direct summa- 
tion over the Matsumara frequencies. In order to re- 
duce the memory storage requirement and accelerate 
the computation process, Pao and Bickers developed 
a renormalization-group computation method 3]. The 
method is based on the assumption that the Green's func- 
tion depends approximately on T only through the Mat- 
subara frequency. The computation starts at high tem- 
perature To to solve the equations of the Green's function 
at selected numbers {n} with cutoff Ao for the Matsubara 
frequencies. Since the Green's function decreases with 
Matsubara frequency very fast at high temperature, the 
selected numbers is not necessarily too many. Then at 
lower temperature T±, the selected numbers correspond 
to lower frequencies. The equations for the functions at 
these lower frequencies are solved and the functions at 
some high frequencies ui n > ojn are approximated as the 
ones calculated at To. For example, G{k,iLo n )T 1 ,n>N a ~ 
G(k,iu' n ) Tg with ui n = TrT 1 (2n+l) = nT (2n' + \) where 
n' is the selected number. The series summation in Eq. 
([I]) is carried out using the staircase rule. That is, sum- 
mation of f(n) over the range n\ < n < ri2 — 1 with ri\ 
and ri2 the any two nearest-neighbor selected numbers 
is given by f (n\){ni — n\). By repeating this sequence, 
the equations for determining the Green's function are so 
solved at lower temperatures. 

The key problem in solving such integral equations by 
the direct summation treatment is how to accurately take 
the series summation with the selected numbers. For 
numerically computing the series summation, 



(3) 



n=0 



S(r,r) = ~v eff (r,T)G(r,T). 



(2) 



the present author has introduced an algorithm that 
sums over selected numbers with corresponding weights. 
The basic idea of this method is described as following. 
Suppose / as a function of continuously variable x is lo- 
cally smooth. Between the selected numbers n\ and 



2 



with n 3 = ni + 2h and h an integer, f(n) can be expanded 
as 



/(n) w /(ni) + ci(n - ni) + c 2 (n - ni) 2 



(4) 



where the coefficients Ci and c 2 are determined by the 
function values /(n 2 ) with 77 2 = ni + ft (the midpoint 
between ri\ and 713 also selected) and f(n 3 ). They are 
given by 



ci = [-3/(m) + 4/(n 2 ) - f{n 3 )]/2h, 
c 2 = [f(n 1 )~2f(n 2 ) + f(n 3 )]/2h 2 . 



Then using the results, 



J2 k = n(n + l)/2, 



^k 2 = n(n+l)(2n+l)/6, 



(5) 
(G) 



(7) 
(8) 



the summation of f(n) over the range n\ < n < n 3 — 1 
is obtained approximately in terms of /(tii), f(n 2 ), and 
fins) the values of / all at the selected points; the coef- 
ficients attached respectively to these values are the cor- 
responding weights depending only on h. By repeating 
this procedure to a large cutoff number, the summation 
in Eq. ((3]) is then obtained. The algorithm is proved to 
be very accurate. 

Though the convolution with the discrete numbers can 
be treated as series summation, numerical computation 
with the discrete Fourier transform is much easier. It is 
even faster provided the transform is performed using a 
high efficiency algorithm. In this work, we will develop 
an algorithm to the discrete Fourier transform. The accu- 
racy and efficiency of the new algorithm will be justified 
with examples. 

In the later part of this paper, we will apply 
the algorithm to the physical problem studying three- 
dimensional interacting electron gas (3DEG) under the 
renormalized-ring-diagram approximation (RRDA) [f| 
and compare the ground-state energy so obtained with 
existing results of the Monte Carlo (MC) simulation [f| 
and the random-phase approximation (RPA). RRDA sat- 
isfies the microscopic conservation laws [H, Q • It has not 
so far been applied to 3DEG because of the numerical dif- 
ficulty in self-consistcntly solving the integral equations 
determining the Green's function. 



II. APPROXIMATION FOR THE DISCRETE 
FOURIER TRANSFORM 

We here consider the discrete Fourier transform 

rib 

F(k)= E /(?')exp(-ifei) (9) 

j=n a 

where f(j) is defined in the range n a < j < rib with n a 
and rib being integer numbers and A: is a real parameter 



in the range (— 7r, 7r). To find out an approximation for 
it, we firstly analyze the following summation in small 
range (711,713) with 773 — 771 = 2h and ni, 713 and h all 
integers, 



713—1 



F(n 1 ,n 3 ;k) = ^ f(j) exp(-ikj). 



(10) 



For large k, since exp(— ikj) is a rapid oscillating factor, 
f(j) exp(-ikj) cannot be regarded as a smooth function 
of j and the previous algorithm cannot be applied here. 
However, for smooth function f{x) in the range ri\ < x < 
773 , f{j) can be expanded as in Eq. We can then 

obtain an approximated result for F{rii, 77,3; k). We need 
the following summation 



713 — 1 



3=ni 



e~xp(—ikrii) 



1 — exp(— i2kh) 



1 — exp(— ik) 
= exp(— ikni)y(k) (11) 

with y(k) = [1 — exp(— i2kh)]/[\ — exp(— ik)]. Then we 
have 

713 — 1 

S2(k) = ^2 ( j - 77i)exp(-7fcj) 

j=ni 

= iexp(—ikni)dy(k)/dk (12) 

"3- 1 

Sa(k) = C? _ m) 2 exp{~ikj) 

= -expi-ikn^cfyity/dk 2 . (13) 
Substituting Eqs. flU and (TJJl-CCI]) into Eq. ^U\), we get 
F(m,n 3 ;k) ^ fin^Sxik) + Cl S 2 (k) + c 2 S 3 (k). (14) 
Using Eqs. (J5J) and ([B]), we obtain 

F(t7i,77 3 ; fc) W 7«l(fc)/(71l) exp(-ifc77i) 

+W2(k)f(n 2 ) cxp(~ikn 2 ) 
+w 3 (k)f(n 3 )exp(-ikn 3 ) (15) 

where the weight functions wi_2, 3 (k) are given by 

3 dy(k) 1 d 2 y{k) 



wi(k) = y(k) - i 



w 2 {k) = exj>(ikh)[i 
w 3 (k) — exp(i2kh)[ 



2h dk 2h 2 dk 2 ' 
2dy(k) 1 d 2 y(k) 



h dk h 2 dk 2 " 
1 dy(k) 1 d 2 y{k) 
2h dk 2h 2 dk 2 1 



Note that these weights depend only on the parameters 
k and h. Equation (|15p means that the summation over 
the range 7i! < j ' < n 3 — 1 can be approximately obtained 
by the three values of the function at 771, 712 and 773. 



3 



Now, we select equal spaced numbers 
(tii, 7i2j • • • ) ri 2m+i) with integer stride h and con- 
sider the summation 



»2m+l-l 



F(n 1 ,n 2m +i;k) = ^ f (j) exp(-ikj) . (16) 
Using Eq. (|T5|) . it can be expressed as 

m 



—ikri2i 



[Wi 



n 2 i-\)e 



-ikn 2 e-i 



1=1 



+w 3 (k)[f(n 2m+1 )e-* kn ^ - f( ni )e~* kn i} 
= w e (k)F e (k) + w {k)F {k) + w x (k)F x (k) (17) 



with 



Fz(fc) = /(n 2m+1 )e- lfc ™ 2 ™+ 1 - /Me-*" 1 (18) 
Fo(k) = ^/("M-Oe - **""- 1 +^x(fc) (19) 



^e(fc) = ^/(n 2 ,)e" 



(20) 



and w e (k) = w 2 {k), w (k) — wx(k) + Wa(k) and w x (fc) = 
[w 3 (fc) - wi(fc)]/2 = Wi(fc) - 1/2. Clearly, F D (k) is the 
summation over the odd terms with the trapezoid rule 
which counts half the values at the two ends, F e (k) is 
that over the even terms, and F x (k) is the extra term. 
Inserting the function y(k) in the formula Wi r 2,3(k), the 
final expressions for w e (k), w (k), and Wi(k) are obtained 
as 

w ^ k ) = [ , bm ^L -2cos(fc/i)]//t(l - cosfc), (21) 



w (k) = 
Wi{k) = 



L /itan(fc/2) 

sin 2 (kh) 

— — - w e [k) cos(kh), 

n{± — cos k) 

h sin k — sin(kh) cos(kh) 1 



(22) 



2Ml-COsfc) ^e(k)Mkh). 

(23) 



At k = 0, these functions take their limit values 

MO) = |(2+ -^), 
u»i(0) = 0, 

which are the same weights obtained previously for the 
series summation 0. The formulas ([!?]) -([2"5]) are the 



main result here. Equation (fT7|) is valid when the func- 
tion f{x) is smooth in each segment (n 2 ^_i, n 2 £ + i). As 
shown above, only in each of these segments, f(x) is ap- 
proximated as a parabolic function and then the sum- 
mation over the terms with the oscillating exponential 
factor is performed exactly. By separating the whole 
range (n a ,nb) into several pieces with different strides 
according to the behavior of the function and using the 
summation rule given above in each piece, we then ob- 
tain the discrete Fourier transform. This reduces greatly 
the memory storage and enhances the computation ef- 
ficiency. Here is a remark: Except for the number of 
the selected points is odd in each piece, there is no con- 
straint on the numbers of values of the input function / 
and the output results and no constraint on the relation 
between the stride of k and the total number of the se- 
lected points {fii}, which is different from the condition 
of the fast Fourier transform. Therefore, it is convenient 
for using. 

Some problems may be related to the sine or cosine 
transform: 



»2m+l- 1 



j=ni 

Il2m+1- 1 

C(k)= fU)™s(kj). 



(24) 



(25) 



By recognizing the functions cos(fcj) and sin(fcj) are 
respectively the real and negative imaginary parts of 
e~xp(—ikj), from Eq. (I17|) . we obtain 

S(k) = w e (k)S e (k) + w (k)S (k) - is s (fc) - Wi (k)C x (k) 
C(k) = w e (k)C e (k)+w (k)C (k)-^C x (k) + Wl (k)S x {k) 



with 



S x (k) = /(n 2m+ i)sin(fcn 2m+ i) - /(n 1 )sin(fcn 1 ) 

rn ^ 

S o{k) = ^ f(n 2 e-i) sin(fcn 2 i-i) + -S x (k) 



m 

S e [k) = ^ f(n 2 e) sm(kn 2e ) 

C x {k) = /(n 2 m+i)cos(fcn 2m +i) - /(ni)cos(fcni) 

rn _^ 

c o(k) = /(n 2 i-i) cos(fcn 2 i-i) + ^C x (k) 

I) i ~ 



rn 



c e{k) = y^/(n 2 <)cos(A-n 2 /). 
t=i 

To test the accuracy and efficiency of the summation 
rule, we here consider an example, 

oo 

U ( x ) = V 2 + L n)2 expH2™a:) (26) 
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FIG. 1: (color online) H(x) as function of x at parameters 
p = 5 and p = 10°. Solid lines represent the exact function. 
Circles are the numerical results obtained using the present 
summation rule. 




FIG. 2: (color online) g(r) as function of r (normalized by 
P = 1/T) at parameters T/£ = 0.01 0.1 and 0.3. Solid lines 
represent the exact function. Symbols are the numerical re- 
sults obtained using the present summation rule. 



where p is a parameter. The exact result of the summa- 
tion is 

1 exp(-px) + cxp[p(x - 1)] 

n W = o 1 1 — ^ ' 

2 1 — cxp(— p) 

for < x < 1. II(x) is a periodic function of x with pe- 
riodicity 1. For numerical calculation, because the term 
under summation is even with n — > — n, we rewrite Eq. 
(1261) as 



1 OO 

n(x) = — + 2j2 

P to 



P 



P 



(2nn) 



■cos(2nnx). (27) 



To get the summation converged in numerical calcula- 
tion, the cutoff Nq that is the terms we need to sum 
should be much larger than p/2ir. For example, suppos- 
ing p = 10 5 , we may take N = 100p/2vr « 1.6 X 10 6 . 
For x ~ 0, the contribution from the remaining term 
of n > N is 0(p/2n 2 No) ~ O(1/100tt). According to 
the summation rule, instead of summing term by term 
within the cutoff, summation in Eq. (|27[) is taken over 
only selected numbers. We here use our previous number- 
selection scheme p|: The selected numbers distribute 
in L successively connected blocks (pieces) in the pos- 
itive integer-number axis. Each block contains M equal 
spaced numbers (selected). (Including the two ends, 
there are M + 1 numbers in each block. The number M 
is redefined here. In different from the previous notation 
where M was defined as the total number including the 
two ends, here it counts the number on one end.) The 
stride (or the length between two selected numbers) in 
the £th block is he — h 1 with h a constant integer num- 
ber. For this example, we here use [h,L,M] = [2, 19,4]. 
The total number of the selected numbers is LM+1 = 77, 
but the cutoff is N = Mh L /(h-l)-(M-h+l)/(h-l) « 
Mh L /(h - 1) = 2 21 w 2.1 x 10 6 . Therefore, the cutoff 
should be large enough for p < 10 5 . Figure 1 shows the 



present numerical results (circles) for Tl(x) as function 
of x at two parameters p = 5 and p = 10 5 . The calcu- 
lation is compared with the exact formula given by the 
solid lines in Fig. 1. For p = 5, Fig.l just shows the 
result within a periodic range. For p = 10 5 , Tl(x) takes 
sizable value when x is close to or 1 within the peri- 
odic range < x < 1 and is symmetric about x = 0.5. 
Here we plot only the result for x close to for p = 10 5 
because the scale p is too large to depict the result in 
a complete periodic range. It is seen that the numeri- 
cal calculation very accurately reproduces the exact re- 
sults. As for the efficiency, because the summation of 
iVo ~ Mh L /(h — 1) terms is approximately obtained by 
summing over only LM + 1 terms, the efficiency c can be 
defined as N /[LM + 1]. That is, 

c^h L /L(h-l). (28) 

For the present example, the efficiency is c = 2.7 x 10 4 . 

In dealing with a physical problem, we may face to 
summation like 



n(a;,pi,p 2 , 



f(n,Pi,P2, ■■■) exp(-i2irnx) 



where f(n,pi,p2, • • • ) cannot be explicitly expressed but 
given numerically. n(x, p\ , pi , • • • ) as function of its argu- 
ments needs to be computed numerically in the region of 
(x,pi,p2, ■ ■ ■) where the summation converges. Accord- 
ing to the present summation rule, first, we only need to 
numerically calculate / at the selected numbers, which 
saves computation time and memory storage for /. Sec- 
ond, the high efficiency summation algorithm saves time 
significantly for getting Tl(x,pi,p2, • • •)• 

To effectively apply the present summation algorithm, 
we here consider another example, 



9 (r) - i 



OO 

E 



exp(-w n r) 



(29) 



where f3 = 1/T, uj n is the fermionic Matsubara frequency, 
and T,(iui n ) = A 2 /(iu n + £) is the self-energy with A a 
parameter. The exact result is 



l(r) 



\(1 + ^)F(-E)exp(-Er) 



for < t < /3, where E = v^ 2 + A 2 and F(E) = 
l/[exp(/3_E) + 1] is the Fermi distribution function. Note 
that G(iu) n ) = l/[iuj n — £ — E(iw„)] ~ l/iu> n as n — >• oo. 
Therefore, the summation in Eq. (f2U]) for r — » is not ab- 
solutely converging but converges conditionally. In such a 
case, one usually makes use of auxiliary function to accel- 
erate the convergence in numerical calculation. We here 
choose the auxiliary function as G°(iuj n ) — l/(iu> n — £). 
The summation 

9°{t) = - ^2 G°(iuj n )exp(~iuj n T) 

n — — oo 

= -F(-S)exp(-£r), for0<r</3 

is known. We then need to do the numerical calculation 
given as, 

5( T ) = g E *e{[G(iu, n ) - ^(ioj^e-^} + 9 °(t) 

where the summation is absolutely converging in the 
limit t — y 0. Using our number-selection scheme with 
[h,L,M] = [2, 17,4], we numerically calculate g(r). The 
obtained results (symbols) are shown in Fig. 2 for pa- 
rameters T/i = 0.01, 0.1, and 0.3 and A/£ = 0.5. The 
solid lines in Fig. 2 represent the exact formula of g(r). 
Clearly, the numerical computation is in very good agree- 
ment with the exact result. 



III. 3DEG UNDER RRDA 

We now apply the above algorithm to study 3DEG un- 
der RRDA. The uniform 3DEG is a fundamental system 
in solid-state physics [7[ and has been extensively stud- 
ied for developing the exchange-correlation functional of 
local density in the frame work of density-functional the- 
ory Q. Within the Green's function approach, most of 
the existing works for studyin g th e system are based 
on perturbation expansions |9l-fl2l]. RRDA is consid- 
ered to be superior to perturbation expansions because 
it satisfies the microscopic conservation laws 1, 2\. The 
similar scheme, the fluctuation-exchange approximation, 
has been extensively applied to the Hubbard models 
for studying the mechanism of high-temperature super- 
conductivity in cuprates [1, Il3l - [20l ]. In different from 
the Hubbard models that describe narrow-band electrons 
with short-range Coulomb repulsion, the 3DEG is a sys- 
tem of electrons with infinitive band width and long- 
range Coulomb interactions. The energy scale of an elec- 
tron in 3DEG is much larger than that in the Hubbard 



O = 



f \aaa 



FIG. 3: (color online) Diagrammatic expressions for 'free en- 
ergy' functional <3>. The line represents the Green's function 
and the wavy line is the Coulomb interaction. 



model and one has to treat the summation over Matsub- 
ara frequencies with a much larger cutoff in the numer- 
ical computation. Since the Green's function needs to 
be self-consistently determined by coupled integral equa- 
tions and the numerical computation is not easy without 
special method, RRDA has not been applied to 3DEG. 
Our objects here are to test the efficiency of the present 
numerical algorithm and to examine the applicability of 
RRDA to 3DEG. 

The three-dimensional electron system with density n 
at temperature T is embedded in a uniform neutralizing 
background of positive charge. The Hamiltonian of the 
system is given by 



H 



Ckcr 



1 

2V 



E v ^ c l 



k+qcr k' 



-qa 



Ck'cr'Cko 



kk' qcrcr' 



where Ck a annihilates an electron of momentum k and 
spin a, e(k) = fc 2 /2 is the kinetic energy, v(q) = Aire 2 /q 2 
is the Coulomb interaction, V is the volume of the sys- 
tem, and the term of q = is excluded from the summa- 
tion because of the neutralizing background. Through 
out the paper, we will use the units in which h = ks — 
m = a = e = 1 with m the mass of the electron and a the 
Wigner-Seitz radius. The Coulomb coupling strength is 
characterized by the parameter 



a/as 



(30) 



with as the Bohr radius. The Fermi degeneracy is mea- 
sured by the ratio T/Ep with Ep = kp/2 the Fermi 
energy and kF = (S^n) 1 / 3, the Fermi wavenumber. 

According to the quantum many-body theory, we start 
with the electronic Green's function G{k,iuj n ). It is re- 
lated to the self-energy E(fc, icj n ) via 



G(k,iui n ) = [iui n - £ fc - T,(k,iu n )]~ 



(31) 



where = e(k) — fi and fi is the chemical potential de- 
termined by 



2T 



n - 



— ^2 iuJ n) exp(iuj n rj) 



(32) 



kn 



with rj an infinitesimal small positive constant. With a 
conserving approximation, the self-energy E is given as 
the functional derivative 1211 



E = 6&/6G 



(33) 



6 



where $ is the 'free energy' functional of the system. The 
Green's function so determined satisfies the microscopic 
conservation laws. Under RRDA, $ is diagrammatically 
given by Fig. 3. In space of (r, r) with < t < j3, £ 
reads 



E(r,r) = -G(r,T)W(r,r) 



(34) 



where W(r,r) is an effective interaction between elec- 
trons as mentioned in the section of Introduction. In 
space of (q,iis m ) with u m = 2mirT the bosonic Matsub- 
ara frequency, W is expressed as 



W(q, iv m ) 



v(q) 



1 - v(q)x{q,iVn 



(35) 



where x(<L iv m ) is the bubble as shown in Fig. 3. In 
terms of G, \ i n ( r ? T)-space is given by 



X(r,r) = 2G(r,T)G(-r,- T ) 

= -2G(r,r)G(r,(3-r) 



(36) 



where the use of G(— r, — r) = G(r, — t) = —G(r,/3 — 
t) has been made. The Green's function G is self- 
consistently determined by Eqs. ([3"Tj) and (j3"2")) and (|3"4")l - 
(1551) . These integral equations are solved by iterations. 
Clearly, in each iteration, we need to Fourier transform 
G(fc, iui n ) to G(r, r) to calculate x( r i T )i transform x(r, r) 
to x(q, ii/ m ) to obtain W(q, iv m ), transform W(q, iv m ) to 
W(r,r) to get E(r, r), and finally transform S(r, r) to 
E(fc, iw n ) to return to G(k,iu n ). 

To numerically do the transforms with guaranteed ac- 
curacy, we need to use auxiliary functions. The key 
points about the transforms and the auxiliary functions 
are illustrated below. 

(i) For transforming G(k,iuj n ) to G(r, r), we choose 
G (k,iu) n ) = l/(iw n - £°), with £g = e(fc) - /i and (j,q 
the chemical potential of the non- interacting electron gas, 
as the auxiliary function as did in the second example 
given by Eq. (JMI- G°(*,r) = exp(-^r) for 

< r < f3 has been given in the example. From G°(k, t) 
to G° (r, r) , we need to carefully carry out the integral 



G°(r,r) = 



1 



2tt 2 



dkkF(-e k )eM-&) sin(fcr). 

(37) 



Since it cannot be integrated out analytically, one has 
to integrate it out numerically. Note that for small r ~ 
0, the factor kF(—^) exp(— £°t) decays slowly as k — )• 
oo. In this case, integrating out the slow decaying part 
analytically, we have 



G°(r,r) = 



f 



27T 2 7 



dkkF(e k )exp(-e k T)sm(kr) 



exp(^ r - r 2 /2r) 
(2ttt) 3 / 2 



(38) 



and the remaining integral is performed numerically with 
the Filon's rule [2l| . This formula is useful only for small 




FIG. 4: (color online) Function rG(r, r) at T/E F 
r s = 5. 



0.05 and 



r. For large r, the factor exp(^oT) is large and the round- 
ing error is big. There is the same factor in the inte- 
grand in Eq. (f3"5]). The formula should be reformed as a 
summation multiplied by this common factor. However, 
Eq. (|37p is a proper formula for numerical integration for 
large r. Having G°(r, r) so obtained, we then need to nu- 
merically transform SG(k,iuj n ) = G(fc,iw„) — G°{k,iu> n ) 
to SG(r,r). We select the fermionic Matsubara frequen- 
cies using the parameters [h,L,M] = [2,17,8]. The 
obtained results are almost the same as that using the 
parameters [h,L,M] = [2,17,12]. For the momentum 
k with a cutoff k c — 25/a, its range is separated into 
four regions: (0,kp — A) with A = mm(2T/kp, fc^/3), 
(k F - A, k F + A), (k F + A, 10/a) and (10/a, k c ) and we 
use 100 uniform meshes in each of them. Since the Fermi 
distribution is sharp at k F at low temperature, we need 
to put dense points there. The r range is divided into 
three regions (0,a), (a, 10a) and (10a, 40a) with respec- 
tively 50, 100, and 150 equal-mesh grids. For the range of 
< r < (3, it is divided into 24 segments symmetrically 
about j3/2. The boundary points of the segments in left 
of /3/2 are given by 



T12 = /3/2, 

Tj = Tj + l/4, 

T = 0. 



for j = 1, 
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In each segment (jj,Tj + i), there are 20 uniform meshes 
for coordinating the Green's function G. It is seen that 
the meshes become dense as r — » or — > (3. This choice is 
necessary at low temperature because the Green's func- 
tion G(fc, t) varies dramatically as r — 5- 0. 

(ii) With G°, the natural auxiliary function for trans- 
forming x{ r i T ) to xilii^m) is chosen as x°( r i T ) = 
-2G°(r,r)G°(r, i 9 - t), and x°{li iv m) is given by 



1 



'9 Jo 



dk^p-J(k,q,v m ) (39) 



7 





FIG. 5: (color online) Real part self-energy E r (fc, ioj n ) at 
T/E F = 0.05 and r s = 5. 



with 



J{k,q, v) 



2 (fc+fc - +l//9)lQ fe^ 



7? 2 



— fcg + ^(arctan 



gfc_ 



arctan ) 

v 



and A;± = fc ± g/2. The remaining transform 5%(r, t) = 
%(r, r) — x°( r 7 T ) to #x(°> *^m) is carried out numerically 
using the Filon's rule. 

(hi) From W(q,iv m ) to W(r,r), we need numerically 
transform 8W{q,iv m ) = W(q,iv m ) — v {q) to SW(r,r) 
while v(q) to v(r) is trivial. At large q or large v m , 
5W(q,iVm) becomes v 2 (q)x(q,iv m )- To image the be- 
havior of x{l^ v rn) m these limits, we consider 



2-K<?{q) + v% 



for \e(q) + iv m \ 



— s> oo 



as a measure of it. The transform of e(g)/[e 2 (g) + isf n ] 
from v m space to r space is the same as in the first ex- 
ample given by Eq. ([26| by noting p = e(q)/T. For large 
q and low T, e{q)/T can be very large. Therefore, the 
cutoff for v m should be large enough. We use the pa- 
rameters [h,L,M] — [2,22,8] for selecting z/ m 's (giving 
rise to almost the same results as that of M = 12). The 
cutoff v m is v c = 2 25 irT ss 3.3 x 10 7 7rT. In our calcula- 
tion, we first transform 5W{q, iv m ) to 8W(r, w m ) choos- 
ing SW°(q,iv m ) = v 2 (q)x{0 7 iv m )/[l - v(q)x(0,iv m )] as 
the auxiliary function. 6W (r,iv m ) is given by 

6W°(r,iv m ) = [exp(-q m r) - l]/r 



with q rn = y/—4irx(Qi i^mj- The q integral in the nu- 
merical transform is performed with the Filon's rule us- 
ing 100, 200, and 100 uniform meshes in ranges (0, 1/a), 
(1/a, 10/a) and (10/a, 35/a), respectively. Finally, we 
transform 6W(r,iv m ) to SW(r, r). 
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FIG. 6: (color online) Imaginary part self-energy E r (fc, iw n ) 
at T/Ef = 0.05 and r s = 5. 



(iv) The self-energy is separated to the Fock (that is 
independent of w„) and the remaining terms. The two 
terms are transformed separately from (r, r) space to 
(k, ioj n ) space. 

Note here that under RRDA x(0, iv m ) is not zero 
and is different from RPA. To see it, by inserting the 
RPA self-energy into the Green's function and calculating 
xilji^m), one then finds x(0,ii/ m ) ^ 0. RRDA is such 
a process that the Green's function is corrected again 
and again until the self-consistency is satisfied. A re- 
lated problem is the plasmon excitation in the system. 
With RPA, the frequency of plasmon is determined by 
the singularity of the summation of ring diagrams. Under 
RRDA, however, it should be determined by the singu- 
larity of a two particle propagator. The effective particle- 
hole interaction in the two-particle propagator is deter- 
mined by second functional derivative of $ with respect 
to the single-particle Green's function G. Therefore, the 
diagram of the irreducible two-particle propagator is not 
a simple bubble. 

The maximum value of e(g) here takes a role of cri- 
terion in determining the cutoff v c of Matsubara fre- 
quency v m . The largest q is 35/a, leading to largest 
e(q) = 0.27 x 35 2 E F « 333£ F - For T/E F = 0.01, we 
have v c ~ 10 6 Ep » e(q). In a narrow-band system such 
as Hubbard model, instead of e(g) as appeared above, 
we may use the band width to estimate the cutoff. To 
see this, we start from the more general expression for 
Xo(q,iVm) 



4 A(k,g)[F(e k )-F(e k+q )] 

vl + A*(k,q) 



with A(k,q) — e(k) — e(fc + q). Clearly, v c should be 
much larger than the maximum value of A(fc, q). In the 
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FIG. 7: (color online) Chemical potential fi as function of 
temperature T for parameters, from top, r a =1, 2, 3, 4, 5, 6, 
8, 10, 15, 20, 25, 30. The inset is a zoom in fi — T picture at 
r s = 5. 



FIG. 8: (color online) Inverse compressibility k 1 as function 
of r s at T/Ef = 0.01. The inset shows fx as function r s at 
T/E F = 0.01. 



Hubbard model, the magnitude of the later is in the order 
of the band width. 

The points in the (k, iu) n ) or (q,iv m ) and (r, r) spaces 
may not be necessarily chosen so dense to coordinate the 
functions. Even with so many points, the computation 
time for solving the integral equations of the Green's 
function with a personal microcomputer is only a few 
seconds, which is tolerable to us. 

With the techniques given above, we have solved the 
integral equations. Fig. 4 shows the function rG{r, r) at 
T/Ep = 0.05 and r s = 5. This function varies dramat- 
ically in a region close to origin (0,0) where there is a 
sharp dip in the surface given by rG(r, r). This behavior 
can be seen from part of the free-particle Green's function 
G° (r, t) as given by the last term in Eq. (l38l) . At large 
r, the surface seems like a wave. This wave is related 
to the Friedel oscillations as seen from Eq. (1551) where 
the Fermi distribution function varies drastically at ItF- 
In Figs. 5 and 6, for the same parameters T/Ep = 0.05 
and r s = 5, we show the real part and imaginary part 
of the self-energy T,(k,iu n ), respectively. As u) n — > oo, 
E(fc, iuj n ) goes to the Fock exchange that is real as shown 
in Fig. 5. At ka ~ 2, S r (fc, iw n )ln->oo varies dramatically, 
showing a logarithmic behavior due to the Coulomb in- 
teraction. On the other hand, T,i(k,iuj n )\ n ^ QO becomes 
zero as shown by Fig. 6. In the limit k oo, £(fc,iw n ) 
vanishes. 

Figure 7 shows the chemical potential /x as a function of 
temperature T at various coupling constant r s . At each 
r s , /i seems as a constant at low temperature. Actually, 
/i slightly decreases with T as shown in the inset of Fig. 
7 for r s = 5, which is a general feature of fi because 
electrons occupy higher energy levels with larger density 
of states due to the thermal excitations. 

Shown in Fig. 8 is the inverse compressibility defined 



by 

«-W(|,^-^(gi> T . (40, 

The corresponding /x as a function r s is depicted in the 
inset. At small r s (high density), kT 1 is positive, im- 
plying that the system is stable. While at large r s , k _1 
becomes negative and the system is unstable. The crit- 
ical value is r s w 5 where k goes to infinitive, implying 
that Wigner crystallization may takes place in the sys- 
tem. This critical value r s ss 5 under RRDA drops in 
the range 4.83 < r s < 104 estimated by the earlier works 
[22T - r28| with 4.83 as the prediction of the Hartree-Fock 
perturbation [2j|. The earlier MC result @ for the crit- 
ical value of Wigner crystallization is r s = 67. 

We here need to emphasize that the compressibility k 
is not equal to — x(0, 0)/n 2 . According to the compress- 
ibility sum rule, k is related to an irreducible two-particle 
propagator. As mentioned above, the diagram of the ir- 
reducible two-particle propagator is not a simple bubble 
under RRDA. 

With the results for the Green's function and the self 
energy, we calculate the energy E = (H). E is given by 

E = 5^[2c(fc) + £ a (fc)]n(fc) 

k 

2j^l-v{q)x(s,w m ) [ ' 

where T, x (k) is the Fock exchange part of the self-energy 
^x(k) = S(fc, ioo), and n(k) is the distribution function. 
On the other hand, E can be also obtained from the 
thermodynamic function O as 

E=%-{pn)*v+pN (42) 
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FIG. 9: (color online) Energy Ae and free energy A/ per 
particle as functions of T at r s — 1, 2, 5, and 10. 



where N = nV is the total number of the electrons and 
f2 is given by [3l| 

n = {$ - exp(w„r?)[£G - ln(-G)]}//?. (43) 



Under the conserving approximation, because of Eq. 
(f33|) . the two equations (|4"Tj) and ([42]) are equivalent. 
Having f2, we obtain the free energy F, 



F = Cl + fiN, 



(44) 



which is related to E by F = E — TS with S as the 
entropy. In Fig. 9, we plot the results for Ae = (E — 
E )/N and Af = (F — E )/N as functions of T with 



E /N 



1.105 0.458 



(a.u. 



(45) 



as the ground-state energy in atomic unit (a.u.) given by 
the Hartree-Fock perturbation [29(. E increases with T 
because of the thermal excitations. But the increment is 
less than the heat TS, so F decreases with T. In the limit 
T — > 0, both of them become the same result, the ground 
state energy. Ae and A/ at T — are the ground-state 
correlation energy per particle e c . Shown in Fig. 10 is the 
result for e c as a function of r s . The present calculation 
(RRDA) is compared with MC @ and RPA. By RPA, e c 
is given by [29| 



r RPA 



v(q)x°(Q, il/ )}- 



(46) 



Clearly, in the whole range of the coupling constant plot- 
ted here, the present RRDA calculation is much closer 
to MC than RPA. The present result is even almost the 
same as the MC simulation at r s > 2. Note that the 
result from Eq. (1411) where S^, n(k) and \ are replaced 



FIG. 10: (color online) Ground-state correlation energy e c 
(in atomic unit) per particle as functions of r s . The present 
calculation (circles) is compared with RPA (squares) and MC 
(solid circles). 



with the corresponding functions of free electrons as in 
RPA is not the same as from Eq. (|4"rj|) plus Eq /N because 
RPA is not a conserving approximation. 

In our previous work [5[, we have studied the two- 
dimensional interacting electron gas (2DEG) under 
RRDA. In that time, we could not perform approxi- 
mated discrete Fourier transform from the Matsubara 
frequency axis to the imaginary time axis, but carried 
out the calculation using approximated series summa- 
tion algorithm for the direct summation over the Mat- 
subara frequency. To test the present numerical method 
of the discrete Fourier transform, we have reinvestigated 
the 2DEG system. For solving the integral equations de- 
termining the Green's function, the numerical computa- 
tion with the present algorithm is much easier and faster 
than that with the previous method. Shown in Fig. 11 
is the ground-state energy as a function of the coupling 
constant r s for 2DEG. The present calculation (circles) 
reproduces precisely the previous results (diamonds). For 
comparison, the results of MC simulation (solid circles) 
[3fj| and RPA (squares) are also depicted in Fig. 11. We 
here make a correction for the RPA result in the previous 
work. The previous RPA notation for the ground-state 
energy is not obtained from Eq. (|46l) but from (|4"T|) . That 
is not the usual meaning for the ground-state energy of 
the RPA calculation. 

As seen from Figs. 10 and 11, RRDA reproduces quite 
accurately the ground-state energy of the MC results, 
better for higher dimensional system. A question then 
raises: Why the critical value r s w 5 for the singularity 
of compressibility of 3DEG is very different from the MC 
value 26 for the ferromagnetization or 67 for the Wigner 
crystallization? To answer it, we recall the expression 
/x = (dE/dN)v,T=o and write kT 1 for T = as 



V K dn 2> 



V,T=0- 



(47) 
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FIG. 11: (color online) Ground-state energy e (in atomic 
unit) per particle as functions of r a for two-dimensional elec- 
tron gas. The present calculation (circles) is compared with 
RPA (squares) and MC (solid circles). The diamonds are the 
RRDA results of the previous numerical calculation. 

Clearly, though the ground-state energy is obtained accu- 
rately with an approximation, its second derivative with 
respect to n may in general not be good. Under RRDA, 
the summation of the most divergent ring diagrams gives 
the predominate contribution to the ground-state energy. 
However, to its derivatives, the contribution from the ne- 
glected diagrams may be also significant. 

Having obtained the total energy E and the free energy 
F, we then get the entropy S as 

S=(E-F)/T. (48) 

Shown in Fig. 12 is s = S/N as a function of T at various 
coupling constant r s . The temperature in Fig. 12 is 
normalized with the energy e 2 /a (that is a scale larger 
than Ep for r s > 1.84). As is seen, for smaller r s , s varies 
linearly in a wider range of low temperature. At large r s , 
the tangent of s decreases as T increasing, reflecting the 
strong coupling effect. The low temperature result for 
s at large r s is hard to calculate accurately since the 
difference between E and F is very small. The situation 
has been shown in Fig. 9 where T is normalized with Ep 
that is a smaller scale than e 2 /a for large r s . The overall 
relative error of the numerical calculation for E and F is 
in the range (0.001,0.01). The difference between them 
can fall into the error bars at large r s and at low T. 
The numerical calculation for s is therefore meaningful 
only at high temperature. At large r s , the electrons are 
strongly coupled. The average energy of an electron is in 
the order of e 2 /a. Therefore, e 2 /a is the proper energy 
scale for strong coupling case. 

With the results for entropy S, we calculate the specific 
heat C defined as 

C = T{—) n . (49) 




0.00 0.02 0.04 0.06 0.08 0.10 

T/(e 2 /a) 

FIG. 12: (color online) Entropy s per particle as functions of 
T at r s = 1, 5, 10, 20, and 30. 



In Fig. 13, we plot the result for C as a function of r s 
at various T. C is normalized by C F the specific heat of 
the free electrons, 

C F = tt 2 NT/2E f . (50) 

As seen from Fig. 13, C/C F monotonically decreases 
with r s for r s > 1. For fixed r s , C/C F is smaller at 
higher T reflecting the tangent change of s(T) as shown 
in Fig. 12. The dashed line in Fig. 13 represents the 
Gell-Mann high-density expansion for the specific heat 

m 

c gm /c f = {1 + arsMn/ars) _ 2]/2^}- 1 (51) 

with a = (4/97T) 1 / 3 . C GM jC F first dec reases from unity 
as r s departing from and then turns to increases with 
r s . Of course, it is valid at small r s . As r s — >• 0, the 
result given by RRDA should be close to this expansion. 



IV. CONCLUSION 

By conclusion, we have developed the approximated 
algorithm for the discrete Fourier transform. When the 
function to be transformed is piecewise smooth, its trans- 
formation can be accurately obtained using a number 
of selected points with corresponding weights. This al- 
gorithm reduces the requirement for computer memory 
storage and enhances the numerical computation effi- 
ciency by several orders. Its accuracy has been examined 
by examples. 

We have applied the numerical algorithm to study 
the three-dimensional interacting electron gas under the 
renormalized-ring diagram approximation. The integral 
equations determining the Green's function are easily 
solved by the present numerical algorithm. Since the 
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subara frequency at low temperature in the numerical 
computation is extremely large. With the present algo- 
rithm, instead of computing the functions at about 2 L M 
Matsubara frequencies in each iteration, one needs to cal- 
culate them at only LM + 1 selected ones. The param- 
eters L and M used here are L — 17 (22) for fermions 
(bosons) and M = 8. The requirement for the computer 
memory storage is greatly reduced and the efficiency is 
c = 2 L /L w 7.7 x 10 3 (1.9 x 10 5 ) for L = 17 (22). We 
have obtained the results for the chemical potential, com- 
pressibility, free energy, entropy, and specific heat of the 
system. The ground-state energy obtained by the present 
calculation is in very good agreement with the result of 
Monte Carlo simulation. 



FIG. 13: (color online) Specific heat C normalized by C F that 
of free electrons as functions of r s at T/(e 2 /a) = 0.01, 0.02, 
0.05, and 0.1. The dashed line at small r s is the Gell-Mann 
high-density expansion [32| . 

band width of the system is infinitive, the number corre- 
sponding to the cutoff of the fermionic or bosonic Mat- 
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